#!/bin/bash -e

# based on anno1030.sh, changed filter af from .2 to .15 and .25 to .2.
. $var


echo dbsnp version is: $dbsnp_version


echo java memory: $java_memory

echo $snpeff_db_version


vcf_path=${data_path}/vcf/gatk/

# java -jar snpEff.jar download -v GRCh37.71


# after gatk and snpeff, anno FILTER goes wrong, so run first.
# ${script_path}anno_fp_0424.sh 1.snp.vcf 1.snp.fp.vcf 2>>$log
# ${script_path}anno_fp_0424.sh 1.indel.vcf 1.indel.fp.vcf 2>>$log
	
	
# format_af.sh 1.snp.vcf
# format_af.sh -p2 1.indel.vcf 
# anno_info_vcf.pl 1.snp.vcf 1.af.txt > 1.snp.af.vcf
# anno_info_vcf.pl 1.indel.vcf 2.af.txt > 1.indel.af.vcf

af_anno.sh -p1.snp $1
af_anno.sh -p1.indel $2

echo
echo
echo gatk VariantFiltration
java -Xmx$java_memory -jar $gatk \
	-R $ref_genome \
	-T VariantFiltration \
	-o 1.snp.filtered.vcf \
	--variant 1.snp.af.vcf \
	--filterExpression "FS > 60.0" \
	--filterName "FS_more_than_60" \
	--filterExpression "MQ < 40.0" \
	--filterName "MQ_less_than_40" \
	--filterExpression "HaplotypeScore > 200.0" \
	--filterName "HaplotypeScore_more_than_200" \
	--filterExpression "af < 0.15" \
	--filterName "Allele_frequency_less_than_15%" \
	--filterExpression "DP < 30" \
	--filterName "Depth_less_than_30"
	
echo
echo
echo gatk VariantFiltration
java -Xmx$java_memory -jar $gatk \
	-R $ref_genome \
	-T VariantFiltration \
	-o 1.indel.filtered.vcf \
	--variant 1.indel.af.vcf \
	--filterExpression "FS > 200.0" \
	--filterName "FS_more_than_200" \
	--filterExpression "af < 0.2" \
	--filterName "Allele_frequency_less_than_20%" \
	--filterExpression "DP < 30" \
	--filterName "Depth_less_than_30"



echo
echo
echo gatk CombineVariants
java -Xmx$java_memory -jar $gatk \
	-R $ref_genome \
	-T CombineVariants \
	--variant 1.snp.filtered.vcf \
	--variant 1.indel.filtered.vcf \
	-o 1.snp.indel.filtered.vcf \
	-genotypeMergeOptions UNSORTED

echo
echo
echo snpsift annotate id
time \
java -Xmx$java_memory -jar $snpsift \
annotate \
-id $data_path/ncbi/dbsnp/All_20150605.vcf.gz \
1.snp.indel.filtered.vcf \
> 1.snpsift.dbsnp.vcf \
2>>$log
# -tabix               : VCF database is tabix-indexed. Default: false



echo
echo
echo snpsift annotate id
time \
java -Xmx$java_memory -jar $snpsift \
annotate \
-id ${data_path}/ncbi/clinvar/clinvar_20150106.vcf \
1.snpsift.dbsnp.vcf \
> 1.dbsnp.clinvar.vcf \
2>>$log




echo
echo
echo snpsift annotate cosmic
time \
java -Xmx$java_memory -jar $snpsift \
annotate \
${data_path}/cosmic/CosmicCodingMuts.vcf \
1.dbsnp.clinvar.vcf \
> 1.id.cosmic.vcf \
2>>$log

echo
echo
echo snpsift annotate varType
time \
java -Xmx$java_memory -jar $snpsift \
varType \
1.id.cosmic.vcf \
> 1.cosmic.var_type.vcf \
2>>$log


echo
echo
echo snpsift annotate gwasCat
time \
java -Xmx$java_memory -jar $snpsift \
gwasCat \
-db ${data_path}/snpeff/gwascatalog.txt \
1.cosmic.var_type.vcf \
> 1.var_type.gwas.vcf \
2>>$log



echo
echo
echo snpsift dbnsfp
time \
java -Xmx$java_memory -jar $snpsift \
dbnsfp \
-db ${data_path}/dbnsfp/dbNSFP2.9.txt.gz \
-f Ensembl_transcriptid,Uniprot_acc,Interpro_domain,SIFT_score,Polyphen2_HVAR_score \
1.var_type.gwas.vcf \
> 1.gwas.dbnsfp.vcf \
2>>$log

# echo
# echo
# echo $snpeff ann snpeff2gatk
# time \
# java -Xmx$java_memory -jar $snpeff \
	# -c ${snpeff_path}/snpEff.config \
	# -v \
	# -lof \
	# -noStats \
	# -o gatk \
	# $snpeff_db_version \
	# 1.gwas.dbnsfp.vcf \
	# > 1.dbnsfp.snpeff2gatk.vcf \
	# 2>>$log
	
	# -t 8 \ ??
	
	
	
echo
echo
echo snpeff ann snpeff2gatk
time \
java -Xmx$java_memory -jar $snpeff \
	-c ${snpeff_path}/snpEff.config \
	-v \
	-lof \
	-noStats \
	$snpeff_db_version \
	1.gwas.dbnsfp.vcf \
	> 1.dbnsfp.snpeff2gatk.vcf \
	2>>$log
	
echo
echo
echo VEP annotation no pick '(can not run in background)'
$tools_path/ensembl-tools-release-82/scripts/variant_effect_predictor/variant_effect_predictor.pl \
--offline \
--force_overwrite \
--dir $data_path/ensembl/ \
--fasta $ref_genome \
--refseq \
--hgvs \
--format vcf \
--vcf \
--input_file 1.dbnsfp.snpeff2gatk.vcf \
--output_file 1.vep.vcf




:<<!
echo
echo
echo $gatk VariantFiltration
java -Xmx$java_memory -jar $gatk \
	-R $ref_genome \
	-T VariantFiltration \
	-o 1.gatk.filter.vcf \
	--variant 1.dbnsfp.snpeff2gatk.vcf \
	--filterExpression "QD < 2.0" \
	--filterName "gatk_QD_less_than_2" \
	--filterExpression "MQRankSum < -12.5" \
	--filterName "gatk_MQRankSum_less_than_-12.5" \
	--filterExpression "ReadPosRankSum < -8.0" \
	--filterName "gatk_ReadPosRankSum_less_than_-8" \
	--filterExpression "ReadPosRankSum < -20.0" \
	--filterName "gatk_indel_ReadPosRankSum_less_than_-20" \
	--filterExpression "InbreedingCoeff < -0.8" \
	--filterName "gatk_indel_InbreedingCoeff_less_than_-0.8"
	
!





# cat 1.dbnsfp.snpeff2gatk.vcf \
# | awk -F', ' '/##INFO=<ID=\w+, /{print $1","$2","$3","$4}'

cat 1.vep.vcf \
| awk -F', ' '{if (/##INFO=<ID=\w+, /) { print $1","$2","$3","$4}else{print $0}}' \
> 1.header_modified.vcf



echo
echo
echo snpsift filter
time \
java -Xmx$java_memory -jar $snpsift \
filter \
"(FILTER = 'PASS')" \
1.header_modified.vcf \
> 1.high.vcf \
2>>$log


# echo
# echo
# echo $snpsift filter
# time \
# java -Xmx$java_memory -jar $snpsift \
# filter \
# "(FILTER != 'PASS')" \
# 1.dbnsfp.snpeff2gatk.vcf \
# > 1.low.vcf \
# 2>>$log



$bgzip -c 1.header_modified.vcf > 1.header_modified.vcf.gz
$tabix 1.header_modified.vcf.gz

echo
echo
echo bcftools view for filter FILTER
time \
$bcftools view \
-f "FP_VARIANTCALLER,FP_PRIMERDIMER" \
1.header_modified.vcf.gz \
| bgzip -c \
> 1.fp.vcf.gz \
2>>$log

gunzip -c 1.fp.vcf.gz > 1.fp.vcf

$tabix 1.fp.vcf.gz

$bgzip -c 1.high.vcf > 1.high.vcf.gz
$tabix 1.high.vcf.gz
$bcftools concat -aD 1.high.vcf.gz 1.fp.vcf.gz |bgzip -c > 1.high_fp.vcf.gz
$tabix 1.high_fp.vcf.gz

if [ ! -d isec ]; then
	mkdir isec
fi
$bcftools isec -p ./isec 1.header_modified.vcf.gz 1.high_fp.vcf.gz

cp ./isec/0000.vcf 1.low.vcf

. $cmd_done
